# 	Filename: jurisdictionMI_source.R
# 	Description: Source script to produce results based on regressions where `ratifyrome' is the dependent variable using the multiply imputed data sets and the data set is restricted to aid-eligible states. Called from jurisdiction_driver.R.
# 	Author: Barry Hashimoto
# 	Date: December 2019
# 	For: Barry Hashimoto, "Autocratic Consent to International Law: the Case of the International Criminal Court's Jurisdiction," International Organization.

############################################################################################################

# Set the number of cores for parallel computation. 
NumberCores <- max(availableCores()-2, 1, na.rm = T)

# Package MI data sets. This code is only used for estimating the Firth logits without parallel processing and will not be run otherwise.
if(FALSE){imputed <- mi(full$imputations[[1]],full$imputations[[2]], full$imputations[[3]],full$imputations[[4]],
  full$imputations[[5]],full$imputations[[6]], full$imputations[[7]],full$imputations[[8]],full$imputations[[9]],
  full$imputations[[10]]); all <- imputationList(datasets = imputed)}

# Models of Rome Statute ratification or accession: consent to jurisdiction. Uses mclapply for parallel computing. 
suppressPackageStartupMessages(library(brglm, quietly = T)) # Bias-reduced logit estimation appropriate for cases of separation due to regional and year dummies.

FirthLogit <- function(x){do.call(brglm,args=list(formula=formula,data=x,family=binomial(link="logit"),method="brglm.fit"))}
formula <- as.formula("ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.single + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict + ptssum.mean + medianIMR + democracy*capital.scaled")
model1 <- mclapply(full$imputations, FirthLogit, mc.cores = NumberCores, mc.preschedule = FALSE)
formula <- as.formula("ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.single + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict + ptssum.mean + medianIMR + democracy*demaid.scaled")
model2 <- mclapply(full$imputations, FirthLogit, mc.cores = NumberCores, mc.preschedule = FALSE)
formula <- as.formula("ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.single + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict + ptssum.mean + medianIMR + democracy*multilataid.scaled")
model3 <- mclapply(full$imputations, FirthLogit, mc.cores = NumberCores, mc.preschedule = FALSE)
formula <- as.formula("ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.single + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict + ptssum.mean + medianIMR + democracy*dod.scaled")
model4 <- mclapply(full$imputations, FirthLogit, mc.cores = NumberCores, mc.preschedule = FALSE)
formula <- as.formula("ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.country + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict + ptssum.mean + medianIMR + democracy*capital.scaled")
model11 <- mclapply(full$imputations, FirthLogit, mc.cores = NumberCores, mc.preschedule = FALSE)
formula <- as.formula("ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.country + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict + ptssum.mean + medianIMR + democracy*demaid.scaled")
model22 <- mclapply(full$imputations, FirthLogit, mc.cores = NumberCores, mc.preschedule = FALSE)
formula <- as.formula("ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.country + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict + ptssum.mean + medianIMR + democracy*multilataid.scaled")
model33 <- mclapply(full$imputations, FirthLogit, mc.cores = NumberCores, mc.preschedule = FALSE)
formula <- as.formula("ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.country + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict + ptssum.mean + medianIMR + democracy*dod.scaled")
model44 <- mclapply(full$imputations, FirthLogit, mc.cores = NumberCores, mc.preschedule = FALSE)

# Estimation of the Firth logit models without using mclapply. Not run.
if(FALSE){ 
  model1 <- with(all, brglm(formula = ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.single
  + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict + ptssum.mean + medianIMR
  + democracy*capital.scaled, family = binomial(link = "logit"), method = "brglm.fit"))
  model2 <- with(all, brglm(formula = ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.single
  + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict  + ptssum.mean + medianIMR
  + democracy*demaid.scaled, family = binomial(link = "logit"), method = "brglm.fit"))  
  model3 <- with(all, brglm(formula = ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.single
  + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict + ptssum.mean + medianIMR
  + democracy*multilataid.scaled, family = binomial(link = "logit"), method = "brglm.fit"))
  model4 <- with(all, brglm(formula = ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.single
  + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict + ptssum.mean + medianIMR
  + democracy*dod.scaled, family = binomial(link = "logit"), method = "brglm.fit"))
  model11 <- with(all, brglm(formula = ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.country
  + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict +  ptssum.mean + medianIMR
  + democracy*capital.scaled, family = binomial(link = "logit"), method = "brglm.fit"))
  model22 <- with(all, brglm(formula = ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.country
  + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict  + ptssum.mean + medianIMR
  + democracy*demaid.scaled, family = binomial(link = "logit"), method = "brglm.fit"))
  model33 <- with(all, brglm(formula = ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.country
  + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict + ptssum.mean + medianIMR
  + democracy*multilataid.scaled, family = binomial(link = "logit"), method = "brglm.fit"))  
  model44 <- with(all, brglm(formula = ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.country 
  + common + mixed + islam + gdp.scaled + ruleoflaw + pre98conflict + ptssum.mean + medianIMR
  + democracy*dod.scaled, family = binomial(link = "logit"), method = "brglm.fit"))
  }

############################################################################################################
# Define some functions useful for tabulating the results.

# A function to stop excess verbosity when running matchit().
quiet <- function(x) { 
  sink(tempfile()) 
  on.exit(sink()) 
  invisible(force(x)) 
 } 

# A function to pool the regression results according to Rubin's Rules for multiply imputed data. For use with R models where we get the $coefficients argument in the function's third line, such as glm() and glmer(). Used in the analyzes of unmatched data sets, below.
pool2 <- function(x, n){
  betas <- list()
  for(i in 1:n){betas[[i]] <- summary(x[[i]])$coefficients[,1] }
  betas <- do.call(rbind, betas)
  vmat <- MIextract(x, fun = vcov)
  ses <- list()
  for(i in 1:n){ses[[i]] <-  sqrt(diag(vmat[[i]]))}
  ses <- do.call(rbind, ses)  
  pooled <- mi.meld(betas, ses, byrow = T)	
  pooled <- do.call(rbind, pooled)
  pooled <- t(pooled)
  pooled <- cbind(pooled, dnorm(pooled[,1]/pooled[,2]))
  colnames(pooled) <- c("Est", "SE", "p")
  return(pooled)
	}
	
# A function to pool the regression results according to Rubin's Rules for multiply imputed data. For use with models estimated model "logit" in Zelig. Used in the analyses of matched sets, below.
  pool3 <- function(x, n){
  betas <- list()
  for(i in 1:n){betas[[i]] <- coef(x[[i]])}
  betas <- do.call(rbind, betas)
  ses <- list()
  for(i in 1:n){ses[[i]] <-  sqrt(diag(vcov(x[[i]])[[1]]))}
  ses <- do.call(rbind, ses)  
  pooled <- mi.meld(betas, ses, byrow = T)	
  pooled <- do.call(rbind, pooled)
  pooled <- t(pooled)
  pooled <- cbind(pooled, dnorm(pooled[,1]/pooled[,2]))
  colnames(pooled) <- c("Est", "SE", "p")
  return(pooled)
	}

# A function for determining whether combinations of coefficients are statistically significant using multivariate Normal simulation.
dem.slopes <- function(x,n){
	require(mvtnorm)
	betas <- list()
		for(i in 1:n){betas[[i]] <- summary(x[[i]])$coefficients[,1] }
			variances <- list()
		for(i in 1:n){variances[[i]] <- as.matrix(vcov(x[[i]]))}
	sim.mat <- NULL
		for(j in 1:n){ 
			sim.mat <- rbind(sim.mat, rmvn(100000, betas[[j]], variances[[j]]))
		}
	summed <- rowSums(sim.mat[,(ncol(sim.mat)-1):ncol(sim.mat)])
	estimates <- c(mean(summed), sd(summed), dnorm(mean(summed)/sd(summed)))
	invisible(estimates)
	}

# A function to get the mean and SD of the AIC from logit models.
get.aic <- function(x,n){
	fit <- vector()
	for(i in 1:n){fit[i] <- AIC(logLik(x[[i]]))}
	stats <- c(mean(fit), sd(fit), NA)
	return(stats)
	}
	
# A function to get the mean and standard dev. of the AIC statistics calcualted from logit models run by Zelig.
get.aic.zelig <- function(x,n){
	fit <- vector()
	for(i in 1:n){fit[i] <- from_zelig_model(x[[i]])$aic}
	stats <- c(mean(fit), sd(fit), NA)
	return(stats)
	}

# A function to get model output statistics.
output <- function(x,n){
	mat <- rbind(pool2(x,n)[(nrow(pool2(x,n)) - 1),], 
	               dem.slopes(x,n),
	               get.aic(x,n),
	               c(nrow(full$imputations[[1]]), NA, NA)
	               )
	mat <- round(mat, 3)            
	rownames(mat) <- c("Coef. of Capital in Autocracies", "Coef. of Capital in Democracies", "Model AIC", "N")
	return(mat)
	}

# Create the matrix for the table with tests of hypothesis 1 
col1 <- list(output(model2,10), output(model3,10), output(model4,10), output(model1,10))
col1 <- do.call(rbind, col1)
col2 <- list(output(model22,10), output(model33,10), output(model44,10), output(model11,10))
col2 <- do.call(rbind, col2)
H1table <- cbind(col1, col2)

#print(H1table)

############################################################################################################
# Models with no control variables. Uses mclapply for parallel computing.

FirthLogit <- function(x){do.call(brglm,args=list(formula=formula,data=x,family=binomial(link="logit"),method="brglm.fit"))}
formula <- as.formula("ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.single + democracy*capital.scaled")
model0 <- mclapply(full$imputations, FirthLogit, mc.cores = NumberCores, mc.preschedule = FALSE)
formula <- as.formula("ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.country + democracy*capital.scaled")
model00 <- mclapply(full$imputations, FirthLogit, mc.cores = NumberCores, mc.preschedule = FALSE)


# Method without using mclapply. Not run.

if(FALSE){
model0 <- with(all, brglm(formula = ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.single + democracy*capital.scaled, 
			family = binomial(link = "logit"), method = "brglm.fit"))
model00 <- with(all, brglm(formula = ratifyrome ~ as.factor(region) + as.factor(year) + ratify.trend.country + democracy*capital.scaled, 
			family = binomial(link = "logit"), method = "brglm.fit"))
 } 
 
# Gather model output 
col1 <- list(output(model0,10), output(model00,10))
col1 <- do.call(rbind, col1)

writeLines("")
writeLines("# Print the results for models of whether a leader accepts the jurisdiction of the ICC. These models have fixed region and year effects but NO CONTROL VARIABLES. The first column has models with a shared baseline hazard. First row has the model with a shared baseline hazard. Second row has the model with a country-specific baseline hazard. Main independent variable is Total Aid and Loans.")
print(round(col1, 3))

#suppressPackageStartupMessages(library(memisc, quietly = T))
#writeLines("")
#writeLines("# Reprint the table for LaTex with 3 sig digits.")
#print(toLatex(col1, digits = 3))

# Remove objects that consume lots of memory and aren't needed now.
rm(model0, model1, model2, model3, model4,model00, model11, model22, model33, model44)

############################################################################################################
# Matching with autocracies. Models of Rome Statute ratification or accession: consent to jurisdiction.
  
# Restrict data set to autocracies  
  aut <- full
  for (z in 1:10){aut$imputations[[z]] <- aut$imputations[[z]][aut$imputations[[z]]$democracy==0, ]}

# Discretize capital.scaled at the median.  
  for (z in 1:10){aut$imputations[[z]]$capital.dummy <- ifelse(aut$imputations[[z]]$capital.scaled > quantile(aut$imputations[[z]]$capital.scaled, 0.5), 1, 0)}

# Drop all observations where ratifyrome is NA
  for (z in 1:10){aut$imputations[[z]] <- aut$imputations[[z]][!(is.na(aut$imputations[[z]]$ratifyrome)), ]}

  # Make data and match
  data.list <- list()
  for (i in 1:10){data.list[[i]] <- with(aut$imputations[[i]], data.frame(ratifyrome, gdp.scaled, ruleoflaw, capital.dummy, common, mixed, islam, pre98conflict, ccode, lcode, year, quarter, medianIMR, ptssum.mean, ratify.trend.country, ratify.trend.single))}

  match.list <- list()
  for (i in 1:10){
  set <- aut$imputations[[i]]
  breaks <- 1/2
  gdp.cut <- quantile(set$gdp.scaled, seq(0,1, breaks))
  rol.cut <- quantile(set$ruleoflaw, seq(0,1, breaks))
  imr.cut <- quantile(set$medianIMR, seq(0,1, breaks))
  pts.cut <- quantile(set$ptssum.mean, seq(0,1, breaks))
  match.list[[i]] <- quiet(matchit(capital.dummy ~ gdp.scaled + ruleoflaw + common + mixed + islam 
  + ptssum.mean + medianIMR + pre98conflict, data = data.list[[i]], 
  method = "cem", cutpoints = list(gdp.scaled = gdp.cut, ruleoflaw = rol.cut, ptssum.mean = pts.cut, medianIMR = imr.cut), verbose = F))
  }
  
  aut.match <- match.list # save it
  
# Logit with a common hazard.
  matched.auts <- list()
  for (i in 1:10){
  	matched.auts[[i]] <- zelig(ratifyrome ~ capital.dummy + gdp.scaled + ruleoflaw + ptssum.mean + medianIMR + ratify.trend.single, 
  	model = "logit",  data = match.data(aut.match[[i]]), cite = F)
  }
  
  #pool3(matched.auts,10)
  
# Logit with unit-specific hazards.
 matched.auts.unittrends <- list()
   for (i in 1:10){
   	matched.auts.unittrends[[i]] <- zelig(ratifyrome ~ capital.dummy + gdp.scaled + ruleoflaw +  ptssum.mean + medianIMR + ratify.trend.country, 
   	model = "logit", data = match.data(aut.match[[i]]), cite = F)
   }
   
  #pool3(matched.auts.unittrends,10)


############################################################################################################
# Matching with democracies. Models of Rome Statute ratification or accession: consent to jurisdiction.
  
# Restrict data set to democracies  
  demm <- full
  for (z in 1:10){demm$imputations[[z]] <- demm$imputations[[z]][demm$imputations[[z]]$democracy==1, ]}

# Discretize capital.scaled at the median. 
  for (z in 1:10){demm$imputations[[z]]$capital.dummy <- ifelse(demm$imputations[[z]]$capital.scaled > quantile(demm$imputations[[z]]$capital.scaled, 0.5), 1, 0)}

# Drop all observations where ratifyrome is NA
  for (z in 1:10){demm$imputations[[z]] <- demm$imputations[[z]][!(is.na(demm$imputations[[z]]$ratifyrome)), ]}

# Package data sets
 # imputed.demm <- mi(demm$imputations[[1]],demm$imputations[[2]],
 #   demm$imputations[[3]],demm$imputations[[4]],demm$imputations[[5]],
 #   demm$imputations[[6]],demm$imputations[[7]],demm$imputations[[8]],
 #   demm$imputations[[9]],demm$imputations[[10]])
 # all.demm <- imputationList(datasets = imputed.demm)

# Make data and match. Islam dropped due to low variance.
  data.list <- list()
  for (i in 1:10){
  data.list[[i]] <- with(demm$imputations[[i]], data.frame(ratifyrome, gdp.scaled, ruleoflaw, capital.dummy, common, mixed, islam, pre98conflict, ccode, lcode, year, quarter, quartersinoffice.since1998Q3, medianIMR, ptssum.mean, ratify.trend.single, ratify.trend.country))
  }

  match.list <- list()
  for (i in 1:10){
  set <- demm$imputations[[i]]
  breaks <- 1/2
  gdp.cut <- quantile(set$gdp.scaled, seq(0,1, breaks))
  rol.cut <- quantile(set$ruleoflaw, seq(0,1, breaks))
  imr.cut <- quantile(set$medianIMR, seq(0,1, breaks))
  pts.cut <- quantile(set$ptssum.mean, seq(0,1, breaks))
    match.list[[i]] <- quiet(matchit(capital.dummy ~ gdp.scaled + ruleoflaw + common + mixed + ptssum.mean + medianIMR
    + pre98conflict, data = data.list[[i]], method = "cem", 
    cutpoints = list(gdp.scaled = gdp.cut, ruleoflaw = rol.cut, ptssum.mean = pts.cut, medianIMR = imr.cut), verbose = F))
    }   
    
  demm.match <- match.list

  # Logit with a common hazard. Common and islam dropped due to zero variance between treated and controls.
  matched.dems <- list()
  for (i in 1:10){matched.dems[[i]] <- zelig(ratifyrome ~ capital.dummy + gdp.scaled + ruleoflaw + ptssum.mean + medianIMR + ratify.trend.single, model = "logit", data = match.data(demm.match[[i]]), cite = F)}
  # pool3(matched.dems,10)

  # Logit with unit-specific hazards. All controls except GDP and rule of law are dropped due to model non-convergence when the controls are included.
  matched.dems.unittrends <- list()  
  for (i in 1:10){matched.dems.unittrends[[i]] <- zelig(ratifyrome ~ capital.dummy + gdp.scaled + ruleoflaw + ptssum.mean + medianIMR + ratify.trend.country,model = "logit", data = match.data(demm.match[[i]]), cite = F)}
  # pool3(matched.dems.unittrends,10)
  
######################################################################################################
# See estimates from all models.
# Raw output for models with full data sets
  #pool2(model0, 10)
  #pool2(model1, 10)
  #pool2(model2, 10)
  #pool2(model3, 10)
  #pool2(model4, 10)
  #pool2(model00, 10)
  #pool2(model11, 10)
  #pool2(model22, 10)
  #pool2(model33, 10)
  #pool2(model44, 10)
# Raw output for models with matched data sets
  #pool3(matched.auts,10)
  #pool3(matched.dems,10)
  #pool3(matched.auts.unittrends,10)
  #pool3(matched.dems.unittrends,10)
######################################################################################################
# Package results of regressions using the matched data for tabular presentation.
calcN <- function(x,n){
  vec <- vector()
  for(i in 1:n){vec[i] <- dim(from_zelig_model(x[[i]])$data)[1]}
  return(c(mean(vec), sd(vec), NA))
  }
  
matchresults.col1 <- rbind(pool3(matched.auts,10)[2,1:3], 
							get.aic.zelig(matched.auts,10), 
							calcN(matched.auts,10),
                            pool3(matched.dems,10)[2,1:3], 
                            get.aic.zelig(matched.dems,10), 
                            calcN(matched.dems,10))
                           
matchresults.col2 <- rbind(pool3(matched.auts.unittrends,10)[2,1:3], 
						   get.aic.zelig(matched.auts.unittrends,10),
						   calcN(matched.auts.unittrends,10),
                           pool3(matched.dems.unittrends,5)[2,1:3], 
                           get.aic.zelig(matched.dems.unittrends,5),
                           calcN(matched.dems.unittrends, 10))

matched.results <- cbind(matchresults.col1, matchresults.col2)

rownames(matched.results) <- c("Coef. of Capital in Autocracies", "Model AIC", "N", "Coef. of Capital in Democracies", "Model AIC", "N")

# Make the table of results for regression models of Rome Statute ratification or accession
table <- rbind(H1table, matched.results)

# Round.
table <- round(table, digits = 3)

writeLines("")
writeLines("# Print the table of results for regression models of whether a leader accepts the jurisdiction of the ICC using multiply imputed data.  The first column has models with a shared baseline hazard. The second row has models with country-specific baseline hazards. The first four rows have models with fixed effects & all control variables. Rows 1-4 feature main independent variable is DAC-EC aid, Multilateral aid, Multilateral loans, and Total aid and loans, in that order. Rows five and six use the matched data set using High Capital as the treatment variable. This table setup mirrors that in Table 1 of the main article.")
print(table)

#writeLines("")
#writeLines("# Reprint the table above for LaTex with 3 sig digits.")
#print(toLatex(table, digits = 3))

######################################################################################################
# Balance statistics for autocracies.
 
# Autocracies: Define match.list
 match.list <- aut.match

# Stats for the matched data alone, now
balance.mat.treated <- NULL
for(i in 1:length(match.list)){
balance.mat.treated <- rbind(balance.mat.treated, t(summary(match.list[[i]])$sum.matched[,1]))
}
treated.balance <- cbind(apply(balance.mat.treated,2,mean))
treated.balance <- round(treated.balance, 2)

balance.mat.control <- NULL
for(i in 1:length(match.list)){
balance.mat.control <- rbind(balance.mat.control, t(summary(match.list[[i]])$sum.matched[,2]))
}
control.balance <- cbind(apply(balance.mat.control,2,mean))
control.balance <- round(control.balance, 2)

rownames(treated.balance) <- rownames(control.balance) <- rownames(summary(match.list[[i]])$sum.matched[,1:2])
colnames(treated.balance) <- colnames(control.balance) <- c("Mean")

n.mat <- NULL
for(i in 1:length(match.list)){
n.mat <- rbind(n.mat, summary(match.list[[i]])$nn[2,])
}
N <- c(apply(n.mat,2,mean)[1], apply(n.mat,2,mean)[2])
# Control stats are in the left two columns, then treated
matched.data.stats <- rbind(cbind(control.balance, treated.balance), N)

# Stats for all data, now

balance.mat.treated <- NULL
for(i in 1:length(match.list)){
balance.mat.treated <- rbind(balance.mat.treated, t(summary(match.list[[i]])$sum.all[,1]))
}
treated.balance <- cbind(apply(balance.mat.treated,2,mean))
treated.balance <- round(treated.balance, 2)

balance.mat.control <- NULL
for(i in 1:length(match.list)){
balance.mat.control <- rbind(balance.mat.control, t(summary(match.list[[i]])$sum.all[,2]))
}
control.balance <- cbind(apply(balance.mat.control,2,mean))
control.balance <- round(control.balance, 2)

rownames(treated.balance) <- rownames(control.balance) <- rownames(summary(match.list[[i]])$sum.all[,1:2])
colnames(treated.balance) <- colnames(control.balance) <- c("Mean")

n.mat <- NULL
for(i in 1:length(match.list)){
n.mat <- rbind(n.mat, summary(match.list[[i]])$nn[1,])
}
N <- c(apply(n.mat,2,mean)[1], 	apply(n.mat,2,mean)[2])

# Control stats are in the left two columns, then treated
all.data.stats <- rbind(cbind(control.balance, treated.balance), N)

# Table for matched data
balance.stats.table <- cbind(matched.data.stats, all.data.stats)

writeLines("")
writeLines("# Balance statistics for AUTOCRACIES. Left pair of columns are control and treated for the matched data. Right pair of columns are control and treated for unmatched data.")
print(balance.stats.table)

#writeLines("")
#writeLines("# Reprint the table above for LaTex with 3 sig digits.")
#print(toLatex(balance.stats.table, digits=3))

###########################################################################################################
# Balance statistics for democracies.

# Democracies: Define match.list
match.list <- demm.match

# Stats for the matched data alone, now
balance.mat.treated <- NULL
for(i in 1:length(match.list)){
balance.mat.treated <- rbind(balance.mat.treated, t(summary(match.list[[i]])$sum.matched[,1]))
}
treated.balance <- cbind(apply(balance.mat.treated,2,mean))
treated.balance <- round(treated.balance, 2)

balance.mat.control <- NULL
for(i in 1:length(match.list)){
balance.mat.control <- rbind(balance.mat.control, t(summary(match.list[[i]])$sum.matched[,2]))
}
control.balance <- cbind(apply(balance.mat.control,2,mean))
control.balance <- round(control.balance, 2)

rownames(treated.balance) <- rownames(control.balance) <- rownames(summary(match.list[[i]])$sum.matched[,1:2])
colnames(treated.balance) <- colnames(control.balance) <- c("Mean")

n.mat <- NULL
for(i in 1:length(match.list)){
n.mat <- rbind(n.mat, summary(match.list[[i]])$nn[2,])
}
N <- c(apply(n.mat,2,mean)[1], apply(n.mat,2,mean)[2])
# Control stats are in the left two columns, then treated
matched.data.stats <- rbind(cbind(control.balance, treated.balance), N)

# Stats for all data, now

balance.mat.treated <- NULL
for(i in 1:length(match.list)){
balance.mat.treated <- rbind(balance.mat.treated, t(summary(match.list[[i]])$sum.all[,1]))
}
treated.balance <- cbind(apply(balance.mat.treated,2,mean))
treated.balance <- round(treated.balance, 2)

balance.mat.control <- NULL
for(i in 1:length(match.list)){
balance.mat.control <- rbind(balance.mat.control, t(summary(match.list[[i]])$sum.all[,2]))
}
control.balance <- cbind(apply(balance.mat.control,2,mean))
control.balance <- round(control.balance, 2)

rownames(treated.balance) <- rownames(control.balance) <- rownames(summary(match.list[[i]])$sum.all[,1:2])
colnames(treated.balance) <- colnames(control.balance) <- c("Mean")

n.mat <- NULL
for(i in 1:length(match.list)){
n.mat <- rbind(n.mat, summary(match.list[[i]])$nn[1,])
}
N <- c(apply(n.mat,2,mean)[1], 	apply(n.mat,2,mean)[2])

# Control stats are in the left two columns, then treated
all.data.stats <- rbind(cbind(control.balance, treated.balance), N)

# Matched data stats
balance.stats.table <- cbind(matched.data.stats, all.data.stats)

writeLines("")
writeLines("# Balance statistics for DEMOCRACIES. Left pair of columns are control and treated for the matched data. Right pair of columns are control and treated for unmatched data.")
print(balance.stats.table)

#writeLines("")
#writeLines("# Reprint the table above for LaTex with 3 sig digits.")
#print(toLatex(balance.stats.table, digits=3))


###########################################################################################################
# SATT Statistics for autocracies.
# For details see http://docs.zeligproject.org/articles/getters.html
# and http://docs.zeligproject.org/articles/att.html
 
 unloadNamespace("arm") # This may conflict with sim in the Zelig/Matchit packages.
 suppressPackageStartupMessages(library(dplyr))
 
  REGIME <- 0
  THRESHOLD <- 1/2 # Dichotomize the capital variable at the median.
   
  aut <- full
  
  for (z in 1:10){aut$imputations[[z]]<-aut$imputations[[z]][aut$imputations[[z]]$democracy==REGIME, ]}
  for (z in 1:10){aut$imputations[[z]]$capital.dummy<-ifelse(aut$imputations[[z]]$capital.scaled > quantile(aut$imputations[[z]]$capital.scaled, THRESHOLD), 1, 0)} 
  for (z in 1:10){aut$imputations[[z]] <- aut$imputations[[z]][!(is.na(aut$imputations[[z]]$ratifyrome)),]}
   
  data.list <- list() 
  for (i in 1:10){data.list[[i]]<-with(aut$imputations[[i]],data.frame(ratifyrome,gdp.scaled,ruleoflaw,
  	capital.dummy,common,mixed,islam,pre98conflict,medianIMR,ptssum.mean,ratify.trend.country,ratify.trend.single))}
  	
  match.list <- list()
   for (i in 1:10){
     set <- aut$imputations[[i]]
     breaks <- 1/2
     gdp.cut <- quantile(set$gdp.scaled,seq(0,1,breaks))
     rol.cut <- quantile(set$ruleoflaw,seq(0,1,breaks))
     imr.cut <- quantile(set$medianIMR,seq(0,1,breaks))
     pts.cut <- quantile(set$ptssum.mean,seq(0,1,breaks))
     match.list[[i]]<-quiet(matchit(capital.dummy~gdp.scaled+common+mixed+islam+pre98conflict+ruleoflaw+ptssum.mean+medianIMR,
     data = data.list[[i]],method="cem",cutpoints=list(gdp.scaled=gdp.cut,ruleoflaw=rol.cut,ptssum.mean=pts.cut,medianIMR=imr.cut)))
     }
 
 # Variables where exact balance is achieved between treated and controls are dropped; adding these increases variance but offers zero reduction in bias. See Ho, Imai King, Stuart (2007) and Iacus & King (2011a,b).
 MatchDataList <- lapply(match.list, match.data)
 formula <- as.formula("ratifyrome ~ ratify.trend.country + gdp.scaled + ruleoflaw + ptssum.mean + medianIMR + capital.dummy")
 GetSATT <- function(x){
    do.call(zelig, args = list(formula = formula, model = "logit", data = x, cite = F)) %>%
    ATT(treatment = "capital.dummy", treated = 1, num = 20000) %>%
    get_qi(qi="ATT", xvalue = "TE")
	}

 # Make a vector of all SATT simulations.
 att.pooled <- unlist(mclapply(MatchDataList, GetSATT, mc.cores=NumberCores, mc.set.seed= 100, mc.preschedule=F))

writeLines("")  
writeLines("# Print the 95% interval  plus 25th and 75th percentiles (top) and median absolute deviation (bottom) of the SATT of the high capital variable on consent to the Rome Statute for AUTOCRACIES.")
  print(signif(quantile(att.pooled,c(.025,.25,.5,.75,.975)), digits = 3)) # 95% interval
  print(signif(mad(att.pooled), digits = 3)) # median absolute deviation.
  
writeLines("")
writeLines("# What is the quarterly frequency of ratifications or accessions to the ICC for AUTOCRACIES in the risk set (i.e. autocracies that had not yet ratified)? Print the mean, median, standard deviation, and median absolute deviation in that order.")
average <- data.list[[1]]$ratifyrome # Data is completely observed, hence only one data set used.
	print(signif(mean(average), 2))
	print(signif(median(average), 2))
	print(signif(sd(average), 2))
	print(signif(mad(average), 2))

writeLines("")	
writeLines("# Calculate the above for AUTOCRACIES where capital.dummy = 1 (autocracies with above-the-median capital receipts). Print the mean, median, standard deviation, and median absolute deviation in that order.")
	average <- data.list[[1]]$ratifyrome[data.list[[1]]$capital.dummy==1]
	print(signif(mean(average), 2))
	print(signif(median(average),2))
	print(signif(sd(average), 2))
	print(signif(mad(average), 2))
	

###########################################################################################################
# SATT Statistics for democracies.
 
if(FALSE){   # Code below not run to save time.
 	
  REGIME <- 1
  THRESHOLD <- 1/2 # Dichotomize the capital variable at the median.
    
  demm <- full
  
  for (z in 1:10){demm$imputations[[z]]<-demm$imputations[[z]][demm$imputations[[z]]$democracy==REGIME, ]}
  for (z in 1:10){demm$imputations[[z]]$capital.dummy<-ifelse(demm$imputations[[z]]$capital.scaled > quantile(demm$imputations[[z]]$capital.scaled, THRESHOLD), 1, 0)} 
  for (z in 1:10){demm$imputations[[z]] <- demm$imputations[[z]][!(is.na(demm$imputations[[z]]$ratifyrome)),]}
  
  data.list <- list() 
  for (i in 1:10){data.list[[i]]<-with(demm$imputations[[i]],data.frame(ratifyrome,gdp.scaled,ruleoflaw,
  	capital.dummy,common,mixed,islam,pre98conflict,medianIMR,ptssum.mean,ratify.trend.country,ratify.trend.single))}
  	
  match.list <- list()
   for (i in 1:10){
     set <- demm$imputations[[i]]
     breaks <- 1/2
     gdp.cut <- quantile(set$gdp.scaled,seq(0,1,breaks))
     rol.cut <- quantile(set$ruleoflaw,seq(0,1,breaks))
     imr.cut <- quantile(set$medianIMR,seq(0,1,breaks))
     pts.cut <- quantile(set$ptssum.mean,seq(0,1,breaks))
     match.list[[i]]<-quiet(matchit(capital.dummy~gdp.scaled+common+mixed+islam+pre98conflict+ruleoflaw+ptssum.mean+medianIMR,
     data = data.list[[i]],method="cem",cutpoints=list(gdp.scaled=gdp.cut,ruleoflaw=rol.cut,ptssum.mean=pts.cut,medianIMR=imr.cut)))
     }
     
# Variables where exact balance is achieved between treated and controls are dropped; adding these increases variance but offers zero reduction in bias. See Ho, Imai King, Stuart (2007) and Iacus & King (2011a,b).
 MatchDataList <- lapply(match.list, match.data)
 formula <- as.formula("ratifyrome ~ ratify.trend.country + gdp.scaled + ruleoflaw + ptssum.mean + medianIMR + capital.dummy")
 GetSATT <- function(x){
    do.call(zelig, args = list(formula = formula, model = "logit", data = x, cite = F)) %>%
    ATT(treatment = "capital.dummy", treated = 1, num = 20000) %>%
    get_qi(qi="ATT", xvalue = "TE")
	}

 # Make a vector of all SATT simulations.
 att.pooled <- mclapply(MatchDataList, GetSATT, mc.cores=NumberCores, mc.set.seed= 100, mc.preschedule=F) %>% unlist

writeLines("")  
writeLines("# Print the 95% interval  plus 25th and 75th percentiles (top) and median absolute deviation (bottom) of the SATT of the high capital variable on consent to the Rome Statute for DEMOCRACIES.")
  print(signif(quantile(att.pooled,c(.025,.25,.5,.75,.975)), digits = 2)) # 95% interval
  print(signif(mad(att.pooled), digits = 2)) # median absolute deviation.

writeLines("") 
writeLines("# What is the quarterly frequency of ratifications or accessions to the ICC for DEMOCRACIES in the risk set (i.e. autocracies that had not yet ratified)? Print the mean, median, standard deviation, and median absolute deviation in that order.")
	average <- data.list[[1]]$ratifyrome
	print(signif(mean(average), 2))
	print(signif(median(average),2))
	print(signif(sd(average), 2))
	print(signif(mad(average), 2))

writeLines("")	
writeLines("# Calculate the above for DEMOCRACIES where capital.dummy = 1 (democracies with above-the-median capital receipts). Print the mean, median, standard deviation, and median absolute deviation in that order.")
	average <- data.list[[1]]$ratifyrome[data.list[[1]]$capital.dummy==1]
	print(signif(mean(average), 2))
	print(signif(median(average),2))
	print(signif(sd(average), 2))
	print(signif(mad(average), 2))
	
	}
	
###########################################################################################################
# End of file. Restart R.






